;; Non-robust computational geometry functions
;; Copyright (C) 2021  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.poly
  (:require [mazama.compgeom.core :refer [distance-squared point-orient2d v==]]))

;;;; NOTE
;;;;
;;;; 1. This software may fail if floating point values are used.
;;;;    This is because == is used to compare numbers and no special effort
;;;;    is made to determine if two floating point values are nearly equal.
;;;;
;;;; 2. This software assumes that:
;;;;    the x-axis increases to the right;
;;;;    the y-axis increases DOWNWARD.

(defn length-of-polyline-squared
  "Return the square of the length of the polyline, `S`, a sequence of 2d points."
  [S]
  (if (< (count S) 2)
    0
    (letfn [(f [v [[ax ay] [bx by]]]
              (+ v (distance-squared ax ay bx by)))]
      (reduce f 0 (map vector S (rest S))))))

;;; Reference
;;; https://en.wikipedia.org/wiki/Polygon#Area (but using 1 multiplication
;;; and two additions per term instead of 2 multiplications and 1 subtraction).
(defn area-of-polygon
  "Compute the signed area of the given polygon, `S`, a sequence of 2d points.
  A positive value means the polygon is counterclockwise,
  a negative value means the polygon is clockwise.
  The last point in the polygon can, but need not be, the same as the first.
  If not then the polygon is implicitly closed."
  [S]
  (if (< (count S) 3)
    0
    (letfn [(f [v [[ax ay] [bx by]]] (+ v (* (+ ax bx) (- by ay))))]
      ;; We return the negative of the normal value to account for downward y
      (- (/ (reduce f 0 (map vector S (concat (rest S) [(first S)])))
            2)))))

(defn- polygon==-helper-1
  [S0 S1]
  (every? identity (map v== S0 S1)))

(defn- polygon==-helper-2
  [S0 S1 index]
  (let [S1' (concat (drop index S1) (take index S1))]
    (or (polygon==-helper-1 (rest S0) (rest S1'))
        (polygon==-helper-1 (rest S0) (reverse (rest S1'))))))

(defn- polygon==-helper
  "We assume the polygons, `S0` and `S1`, have the same length and are non-empty."
  [S0 S1]
  (let [item (first S0)
        indices (filter number? (map #(if (v== %1 item) %2 nil) S1 (range)))]
    (some #(polygon==-helper-2 S0 S1 %) indices)))

(defn polygon==
  "Are the two polygons, `S0` and `S1`, ==?
  Regardless of where along the polygons you start and which direction you go."
  [S0 S1]
  (and (== (count S0) (count S1))
       (or (and (empty? S0) (empty? S1))
           (polygon==-helper S0 S1))))

(defn- belt
  "Append the first `n` elements of `coll` to the end of the collection, in order.
  For example, (belt 2 (range 10)) returns (0 1 2 3 4 5 6 7 8 9 0 1)"
  [n coll]
  (concat coll (take n coll)))

(defn convex?
  "Return 1 if `S` is convex and counterclockwise, -1 if convex and clockwise, otherwise false."
  [S]
  (loop [B (belt 2 S)
         status 0]
    (let [[p q r] B]
      (if (nil? r)
        (when (not (== status 0))
          status)
        (let [orient (point-orient2d p q r)]
          (cond
            (== orient 0) (recur (rest B) status)
            (== status 0) (recur (rest B) orient)
            :else (when (== orient status)
                    (recur (rest B) status))))))))
